/*

 Given G(k,t,t') define self energy as a sum of phonon, Coulomb and impurity 
 self energy. All self energy codes reside in this file.

 */

#include <selfenergy.h>

void se::impurity_self_energy(
	const parameters        &params,
	const h5files_container &h5files,
	const int                itau,
	const int                it,
	const arma::vec         &weights,
	complex_matrix          &Vimp,
	complex_matrix_function &G_current,
	complex_matrix_function &impurity_self_energy)
{
	/*
		Formula for the impurity self energy implemented in this function:

		Σ_{nm}(k;t,t') = V_{nn'}(k,k') G_{n'm'}(k';t,t') V_{m'm}(k',k)
		
		This is the second-order Born approximation to the self energy, which is
		second order in the impurity potential. This formula is in terms of the
		very general impurity potential Fourier transform. This potential
		for N impurities in sample area A is given by
		
		V_{nn'}(k,k') = 1/A Sum[ Vi_{nn'}(k,k'), i=1,...,N  ]
		
		The Born self energy then becomes:
		
		
		Σ_{nm}(k;t,t') = N/A^2 V_{nn'}(k,k') G_{n'm'}(k';t,t') V_{m'm}(k',k)
		
		Here N/A = impurity density. Converted to integration over k'

		  Σ_{nm}(k;t,t') 
		= (impurity density)
		 x Integral[ 1/(4pi^2) V_{nn'}(k,k') G_{n'm'}(k';t,t') V_{m'm}(k',k)]

	*/
		
	complex_matrix F;
	int ntmem = params.ntmem;
	int nk    = params.num_mesh_points;
	
	for (int ik=0; ik < params.num_mesh_points; ik++)
	{
		for (int ip=0; ip < params.num_mesh_points; ip++)
		{
			for (int itau=0; itau < params.ntmem; itau++)
			{
				F =   Vimp[ip+ik*nk];
				F = F*G_current[itau + ip*ntmem];
				F = F*Vimp[ik+ip*nk]; 

				impurity_self_energy[itau + ik*ntmem] += F * weights(ip);
				
			}
		}
		
	}

	return;

}

complex_matrix phonon_propagator(double tau, double kmag){
	complex_matrix F;
	F.zeros();
	return F;
}
void se::phonon_self_energy(
	const parameters        &params,
	const h5files_container &h5files,
	const int                itau,
	const int                it,
	const arma::vec         &weights,
	const arma::vec         &tau,
	const arma::vec         &kmag,
	complex_matrix_function &G_current,
	complex_matrix_function &phonon_self_energy)
{


	/* 1. Construct the phonon propagator 

	Σ_{nm}(k;t,t') = D_{nn'}(k,k';t,t') G_{n'm'}(k';t,t') 	
	
	*/

	/* 2. Calculate the Convolution with G_current */
	
	/* 3.  */


	complex_matrix F;
	int ntmem = params.ntmem;
	int nk    = params.num_mesh_points;
	
	for (int ik=0; ik < params.num_mesh_points; ik++)
	{
		for (int ip=0; ip < params.num_mesh_points; ip++)
		{
			for (int itau=0; itau < params.ntmem; itau++)
			{
				F =   phonon_propagator(tau(itau),
				                        kmag(abs(ik-ip)*ntmem));
				F = F*G_current[itau + ip*ntmem];

				phonon_self_energy[itau + ik*ntmem] += F * weights(ip);
				
			}
		}
		
	}
	
	return;

}
